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Abstract. In dense, static, polydisperse granular media under isotropic pressure, the probability density and the correlations 
of particle-wall contact forces are studied. Furthermore, the probability density functions of the populations of pressures 
measured with different sized circular pressure cells is examined. The questions answered are: (i) What is the number of 
contacts that has to be considered so that the measured pressure lies within a certain error margin from its expectation value? 
(ii) What is the statistics of the pressure probability density as function of the size of the pressure cell? Astonishing non-random 
correlations between contact forces are evidenced, which range at least 10 to 15 particle diameter. Finally, an experiment is 
proposed to tackle and better understand this issue. 
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INTRODUCTION 

One of the open issues in the field of disordered, random systems like dense, static granular packings, is the probability 
density of the contact forces and their possible longrange correlations. There is common agreement that the probability 
for large forces decays exponentially H S i S H i, 0, iSSIHllillilHIil 

but the small forces are much 

harder to measure I2l [l6l[l7ll . so that there is still ongoing discussion about the shape of the probability density for 
small forces, possible correlations between the forces, and a predictive model for the force propagation inside dense 
packings of frictional particles ifTsll . 

Furthermore, it is observed that the deformation of particle systems is not affine in general, but displays finite 
distance correlations which are assumed to increase when approaching the jamming transition |6, 19, 20, 21, 22, 2%, 
l24ll25[ I26II27I1 . However, the issue whether these correlations (and possibly anti-correlations due to vortex formation) 
depend on the system size fl^ or not fl^J^d] is not completely resolved yet. 

Since the forces that granular particles exert onto their confinement (walls) strongly fluctuate from one particle to 
the next, so does also the local pressure. When the pressure on the wall is measured with a circular pressure cell, 
performing many independent measurements, one obtains a probability density of the measured pressure, with its first 
moment approaching the mean pressure, and the standard deviation decreasing with increasing cell size. While the 
case of uncoiTelated forces is rather behaving as expected, the case of subtly correlated forces in granular packings 
leads to interesting results and probably can be understood with advanced statistical methods [31, ], which can e.g. 
account for contact number fluctuations. 

Thus, a local measurement of the wall pressure can be far away from the total pressure (or the mean, representative 
pressure), unless "enough" particles are contained in the pressure cell. Besides the question, how much "enough" 
is, also the question of the behavior of the width of the pressure probability density is examined in this study. The 
reason to examine the pressure distribution instead of the force distribution is that the former is much easier to access 
experimentally, as will be outlined at the end of this paper 

In general, more knowledge on the force- and pressure-density functions is needed for the understanding of pressure 
measurements aiming, for example, at a safe design of containers of granular materials such as silos. 



Review on force probability 



The simplest model for the force probability density function is the so-called <7-model, introduced 1995 by Liu, 
et. al. lfl6l l2ll. that describes the occurrence of force chains in disordered geometric packings of granular media on 
the basis of a mean field approximation ifflfisll . In a dense packing of discrete particles, the contact forces that one 
particle in a certain layer feels from above plus its weight, determines the sum of the contact forces on the particles 
below. The magnitude of the two forces at the two contacts below, are the fractions q and 1 — ^ of the sum. In general, 
the weights can also depend on neighboring sites and contacts can open and close. The mean field approximation, 
however, neglects these dependencies and thus simplifies the model vastly. Eq. ([T]i gives us the normalized, scaled 
density function of the inter-particle forces /* =//(/), predicted by the g'-model: 

(^^/f"'exp(-C/*), (1) 

where C is the number of the neighboring particles (below or above). The weakness of the ^-model is the improper 
prediction of the probability to find small forces, see Fig.[T] 




FIGURE 1. Normalized probability density plotted against the normalized force /* = //(/)• The circles are simulation 

data, and the dashed and solid lines are Eqs. (TJ and (|2}, with C = 1.8, respectively. The inset shows the same data in log-scale, 
together with the "quality factor", i.e. the simulation data divided by the fit-function as dots around unity; the fit has deviations less 
than two per-cent for /» < 3 and less than ten per-cent for /» < 5. The g-model, Eq. (Til, is invalid for /, < 0.2 and has varying 
deviations of about ten per-cent for 0.2 < < 5. 



A function, that provides an excellent fit to the simulation data is 



Pf{f*) = 1 -aexp - 



cexp(— t//*) 



(2) 



with the fit-parameters a = 0.983 ± 0.003, b = 0.56 ± 0.05, c = 1 .80 ± 0.02, and d=\OA± 0.7. Note that a similar 
function was found experimentally, see Refs. i32l [33ll . The number of parameters can be reduced by using the 
normalization relations / p fdf = 1 and J f p fdf = I □• The fit-parameters for different densities are given in table [T| 



TABLE 1. Fit results for different volume fractions v, 
where the row e contains the typical relative error of the 
respective fit-parameter. 
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FIGURE 2. The forces of particles that touch the walls fluctuate around a mean value. 

THEORY 

Let us consider a cubical container filled with granular material under isotropic, hydrostatic pressure. Each particle i 
of the Nw particles that touch a wall exerts a force fi (see Fig.|2|i onto the wall and contributes to a (finite: 1 < / < Mv) 
population of forces p{f) with mean jx^y = (/) and standard deviation (7,^, where the subscript w refers to the population 
of all wall-particle forces. If we select one sample j of size n out of this population by applying a circular sensitive 
area D{R) = nR^ that includes the n forces that are acting on this area, we obtain the pressure 

P0)=Pj{n,{R)) ^ ^^^^l^-^y/-'' = (l/«) tfj/Di , (3) 

with the area per particle, D\ = D{R) /n, and n = nj{R). Note that such areas have to be selected such that their center 
- on the selected wall - is at least a distance R away from any other wall. Taking many samples, m':^ \, will result 
in a population of pressure values with probability density pr[P) around the mean jti/> = (Pj{R)) = (1 /m)Y!J=iPi{R), 
with standard deviation 

ap = ^jipj{R))-pLl, 
where the subscript P refers to the population of pressures. 

Central Limit Theorem 

According to the central limit theorem (CLT), and for the corresponding assumptions, the probability density 
function pr{P) of our population of samples Pj{R) provides the same expectation value as for the original population 
/i/> « ^4=1 fi/^^ where A denotes the total surface of the confinement and the sum goes over all particles in the system 



the factor 8 in the denominator leads to some better fit-quality by stretching the Gaussian correction function 



D(R) 



il>(R) 

FIGURE 3. The average force of sample j, divided by the sensitive area D{R) /n per particle, of this detector, leads to a pressure 
Pj{R) corresponding to one measurement j. 



touching the wallsB The CLT also tells us that the probability density function becomes more and more Gaussian: 

2" 
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2\ Op 



(4) 



the larger we chose n, i.e., by increasing the detector size R. The standard deviation then equals Op = Gp{R) = 
{n{R) (7n-) I {D{R)^/n{R)) — a„l {D\ ^/n{R)) =: Pf/y/n, with the pressure that corresponds to the standard deviation 
of the force density function scaled by the area of the pressure cell per particle. 



Confidence Intervals 

Integration from to +°o of the normalized distributions pr{P) and gR{P) gives unity. Now, to gain more 
advanced statistical predictions about the pressure distribution, let us consider a lower and an upper integration limit 
< I^P ™d Zi_a/2 > jJ^p, respectively, such that the integral overEq. dUl equals 

l-a= J gR{P)dP=\- j gR{P)dP- j gR{P)dP=^- 

If the integration limits are chosen such that a fraction a /2 lies outside of the integration range, both to the left and 
the right, the integration limits correspond to the confidence interval 25a — Zi-a/2 ~^a/2^ ^ fraction 1 — o: of the n 
measurements Pj lies within the confidence interval. 

Keeping n — const, and considering a = (probability for finding the measured value in-between our limits then 
equals 1), we expect that the interval of confidence tends to infinity. On the other hand, for a — \ (probability vanishes), 
one expects 5a 0. 

Due to the symmetry of the Gaussian distribution, one can compute 

^a'-= Z\-al2- l^P = l^P-^a/l (6) 
explicitly, using the relation: ierf(l — a) = 5a/(\/2o>) = 5a\/n/ {\/2Pf ) such that: 

« = 2(P^/5„)^ ierf2(l-a) , (7) 

where the pressure Pf = nijCjO^./DiR) = Ok/Di is, of course, a constant for given geometry and p{f). The general 
equation (|7]l gives us information about how many particles n need to be considered in a measuring process, so that a 
measurement Pj{R) lies in-between the error margin ±5^ with probability 1 — a. Note, that 5a and n depend on the 
size R of the sensitive area 5a °^ l/\/n{R) °^ 1 /R. 
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(5) 



^ This is valid under the assumption that all wall-particles are similar, i.e. there is no inhomogeneity in the forces on the wall particles, e.g. as 
function of distance from another wall in the edges of the cuboid. 



Explicit predictions 



q-model 

Now we use the distribution predicted by the q'-model for the analysis. The first two moments obtained from Eq. ([T]l 

1 /2 

are / = (/) and p = ^^{.f)^, which leads to the standard deviation Pq{f), namely aj = (^f^ ~ f^^ = ^1/C(/). 

Considering Eq. ([T]i as a finite population and applying the CLT by taking many samples, one can then replace Pj in 
Eq. dill to get 

= § imRl''''^'^^ " ^ (2/C)ierf2(l - a){Pf/5af . (8) 
Best Fit 

Now we use the distribution as obtained by our fit to the data for the analysis. The first two moments obtained 
from Eq. (|2]i are / « (/) , = 1 .6 1 (/) ^ and the corresponding standard deviation can be computed to aj" —0.78(f). 
In comparison with these results the g'-model provides aj — 0.76 (/) for a dense granular packing with C = 1.8 
corresponding to a special geometry. 

SIMULATION 




FIGURE 4. Snapshot of a dense, polydisperse assembly of N = 8000 particles confined in a cuboid. The volume fraction here is 
V ~ 0.7 and the particle sizes are grey scaled (bright particles are big, dark particles are small). 




FIGURE 5. View on one of the walls for a packing with A' = 20000 particles with v ~ 0.67. (Left) Each circle is a particle in 
contact with the wall; the color coding is the same as in Fig.|4] (Right) Same data as (left), but here each circle radius is scaled with 
the force exerted by the particle on the wall. (Big, bright circles correspond to large forces, whereas small, dark circles correspond 
to small forces on the wall). 

The System 

The systems studied in the following contain perfectly spherical particles with radii r, drawn from a homogeneous 
distribution r,- G ['"min,'"max], where we used rmax/'"min = 2 and 3, and the mass {4-/3)7rrfp, where p represents the 
uniform material density of the paiticles. Since the mass is not relevant in the static limit, we refer to density only with 
the dimensionless volume fraction Vf ~ ^^^4 ''f /I\^pl^i3^ where Lp denotes the length of the simulation volume in 
direction j3 E {x,y,z)- All of the simulations were done in a cuboid volume which is limited by walls that repell touching 
particles. Fig. |4] gives an example of a typical static and dense granular sample, while Fig. |5] shows a representative 
subset of particles that touch a wall (Left) and the corresponding magnitudes of forces (Right), quantified by the radius 
of the circles. 



Molecular Dynamics 

The simulations were performed by means of a molecular dynamics (MD) code in three dimensions, without 
tangential forces like friction. MD simulations are characterized by discretizing time into timesteps At and solving 
the Newton's equations of motion for each particle. In each integration step, the new position of the particles is 
computed from its previous and present positions and accelerations due to forces currently acting on it |34]. MD is 
also referred to as discrete element method (DEM) or as soft-sphere model, i.e. the repulsive forces /" normal to the 
plane through the point of contact depend on the overlaps of the spheres, that replace the contact deformation. The 
linear spring-dashpot (LSD), the non-linear Hertzian model and a hysteretic force model can be used |35], among 
many others. For particle-wall contacts, the same spring constant were used as for contacts between particles. Here, 
only the LSD results are discussed because no qualitative differences could be evidenced for the different contact laws. 



Initial Configurations and Relaxation 

As initial configurations, — 8000 or N = 20000 poly-disperse particles with random initial velocities were placed 
on a regular cubic lattice with low total density. Due to the free space between the particles, the initial order is forgotten 
and the particles disspate energy during collisions. Eventually, the dense, relaxed and disordered granular packing is 
obtained by either applying hydrostatic pressure on the walls or by growing the particles. Due to the dissipative nature 
of the contact law, energy is dissipated and the system reaches a static configuration, where we use as criterion the ratio 



of kinetic and potential energy, £kin/£^pot ~ £, with e < 10^^. If e is small, the particles are typically almost at rest 
and the major contribution to the total energy stems from the contact potential energy between particles and between 
particles and the walls. Since the results were identical for the two preparation procedures, we only mention that there 
are more alternative ways of achieving a static packing, see 1,36,1 . 

RESULTS 

The force- and pressure density functions are gained from the (repulsive) forces between all particles that touch the 
walls. Note that the scaled force-density function for all particles in the bulk was identical to the force-density for 
particle-wall forces, within the considerably larger fluctuations of the latter, only the mean particle-particle force is 
smaller than the mean particle-wall force. Since we are interested in the experimentally accessible pressure measured 
at the wall, we will not present bulk data in the following. 

When a circular measuring area with radius R is put around each wall-particle (for all wall particles with distance 
larger than Rm^x from any other wall), one obtains a set of Pj{R) values and, from these, can compute the mean values 
and standard deviations. A typical set of wall-particles is shown in Fig. |5] Note that Rmax is introduced, so that the 
same set of wall-particles is used for the computation independent of R and, if ^max/'"max becomes larger than 12-13, 
the statistics becomes bad, since too few particles close to the center of the wall are taken into account as the centers 
of pressure cells. 

For each sample with given R, from the corresponding pressures Pj{R), the histograms are obtained, as shown in 
Fig. |6] for three different /^-values. The larger the cell size, the more the pressure density function appears Gaussian 
and, as Fig.|7]shows, the smaller the standard deviation Gp becomes. Note, however, the interesting fact, that the decay 
of o> is steeper than the simple relation CTp 1 /R, expected from the central limit theorem. 
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FIGURE 6. Probability density for the pressures from many samples with sizes R^, = R/rmax = 1-0, 2.5, and 5.0. Points are 
simulation data and lines correspond to Eq. Q The left panel shows the same data as the right, only the latter has a logarithmic 
vertical axis. 

Fig. |7] also contains steps for small cell sizes due to integer jumps in the particle number n{R). For larger R, the 
change of Op becomes smooth and independent of For every kind of distribution of the population we should get 
- according to CLT - a Gaussian distribution of the Pj(R) around the population mean. Eq. (|7]) is a consequence of 
these simplifications. But our simulations show (n) o= 5^''^ (see Fig.[8]l or, what is aequivalent, deviations from the 
Gaussian shape of the curves for large R as well. Note, that for very small radii only a very small amount of particles 
will contribute to the pressure, i.e., at least, the central particle will contribute to the pressure. 

These correlations vanish if we replace the forces in our data files by randomly and uniformly distributed forces. 
These (faked) data then show agreement with Eq. (|7]), see also Fig. [8] as expected, and finally confirm that correlations 
do occur in our systems (samples of particles that contribute to the pressure measured at the walls). 
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FIGURE 8. Number of particles averaged over many samples with given radius R, for many radii drawn from the interval of 
confidence for a = 0, scaled by the population mean. The simulations were performed with A' = 20000 particles and LSD. The data 
with slope of — 1 .5 gained by our MD simulations reveal correlations between the forces of the wall-particles because they do not 
match to theory, see Eq. Q, that predicts a slope of —2. Data with slope of —2 are "faked" by replacing the forces by uniformly 
distributed random forces. 
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Correlation Function 

An alternative way to examine the correlations in a set of forces is to compute the correlation function 

C(r) ^ (9) 

where the fa are the forces of the center particles, for a = i, and of all other particles at distance r in the pressure 
cells, for a — j. As displayed in Fig.|9](Left), the correlations decay about two orders of magnitude within a distance 
of approximately 15 particle (maximal) radii. For small R'^-^^^, this contrasts the results from the confidence interval (or 
standard deviation) measurements, where no change of behavior could be evidenced. These results are not changing 
for larger TJ;;^^. 

The data for different densities in Fig. |9] (Right) show that the correlations decay much faster for larger densities. 
Note that there is no qualitative difference visible when one examines o/jJ-p as function of R* for different densities: 
The absolute values decrease with increasing density, but the slope of —3/2 remains independent of the density (data 
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FIGURE 9. (Right) Correlation function 1 — C{r*) as function of the scaled distance from the center of the pressure cell, 
r* = r/rmax, for different sizes of the pressure cells /fmax' given in the inset. (Right) Correlation function 1 — C(r*) for different 
densities v, as given in the inset. 



not shown). Thus, in conclusion, the pressure cell approach is capable of detecting a different type of correlations, 
which is not caught by the classical correlation function approach. 



CONCLUSIONS 

Our main objective was achieved by determining how the number of particles contributing to the pressure on a given 
sensor correlate with the interval of confidence, which is a measure for the width of the pressure density functions. 
As expected, the confidence interval increases the smaller the detector size is chosen, i.e., the worse the statistics 
becomes. Surprisingly, the relation {n{R)) °^ 5a(/?)"'^ is observed, i.e. for a certain amount of particles, given a 
desired probability 1 — a, the measured pressure values can be expected to be closer to the average, (/?)), as assumed 
from CLT. Thus, our simulations predict a better confidence in measured data. 

As possible reason for this, one has the correlations between forces exerting by close-by the particles on the walls. 
These correlations range over a rather long distance. The rather limited sensor size (and the related fluctuations of 
particle numbers) could excluded as source of this effect, because we found {n{R)) «: 5a{R)~'^ for a fully random 
(uncorrelated), uniform force distribution. 

Finally, we note that our systems were rather small, so that we cannot exclude the possibility that our observations 
are due to a finite size effect. Therefore, much larger simulations should be performed to confirm the present results 
and better understand the source of the correlations as evident from the width of the pressure probability density 
function. Furthermore, we analysed our simulation data only at the walls, so that a direct experimental access to the 
same information is possible. In addition, the pressure correlations presented here should be confronted to bulk stress 
data in order to learn if this is a wall effect or intrinsic also to bulk particulate systems far away from the waUs. 
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